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RESUMEN 

Se describe un algoritmo eficiente tipo Monte Carlo para solucionar una clase restringida de problemas en 
la transferencia de radiation. Esta clase incluye varios problemas de inferos en la astrofisica, incluyendo la 
dispersion de la luz ultravioleta y visible por granos. El algoritmo toma en cuenta la dispersion multiple. Se 
describe el algoritmo, se presentan algunas optimizaciones importantes, y se muestra explicitamente como se 
puede usar el algoritmo para estimar cantidades como la intensidad emergente y promedio. Se presenta dos 
pruebas del algoritmo, se examina la importancia de las optimizaciones, y se muestra que el algoritmo se puede 
aplicar de manera litil a los problemas opticamente delgados, los cuales son a veces considerados limitados a 
aproximaciones explfcitas de dispersion sencilla mas atenuacion. 

ABSTRACT 

We describe an efficient Monte Carlo algorithm for a restricted class of scattering problems in radiation transfer. 
This class includes many astrophysically interesting problems, including the scattering of ultraviolet and visible 
light by grains. The algorithm correctly accounts for multiply-scattered light. We describe the algorithm, 
present a number of important optimizations, and explicity show how the algorithm can be used to estimate 
quantities such as the emergent and mean intensity. We present two test cases, examine the importance of 
the optimizations, and show that this algorithm can be usefully applied to optically-thin problems, a regime 
sometimes considered limited to explicit single-scattering plus attenuation approximations. 
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1. INTRODUCTION 

We present a Monte Carlo algorithm for the numerical solution of a restricted class of scattering problems 
in radiation transfer. This class consists of problems in which the true emissivity (the emissivity of unscattered 
light) and the opacity, albedo, and scattering properties are specified a priori and are zero outside a finite 
domain, but have arbitrary distributions within that finite domain. In particular, no restrictions are placed 
on the scattering phase function or on the symmetry of the problem. Examples of problems that satisfy these 
requirements include those in which the opacity is dominated some combination of grain opacity, electron 
scattering, resonance-line scattering, Rayleigh scattering, or Raman scattering. For example, the transfer of 
ultraviolet and visible light in the presence of grains and of Lya photons in the presence of neutral hydrogen 
and grains. For simplicity, in this paper we do not consider time-dependent problems or polarization, however, 
as we discuss in §|[ the algorithm can be easily extended to cover these cases. 

We have developed this algorithm over the last several years, initially independently (Watson 1994; Henney 
1994ab; Henney 1995; Henney & Axon 1995; Burrows et al. 1996; Sahai et al. 1998; and Henney 1998) 
and then collaboratively (Stapelfeldt et al. 1999; Watson et al. 2001). The algorithm currently incorporates 
several important features and optimizations, including the use of forced scatterings and forced interactions, 
the construction of estimators from unbiased samples of scattering events, and a very efficient integrator for 
optical depth. 

Monte Carlo algorithms have been widely used in radiation transfer (see §||). It seems clear to us that 
aspects of this algorithm have been independently discovered and appear to be common knowledge among 
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TABLE 1 
NOTATION 

Symbol Definition Units 

r Position. cm 

Vi Position of the i-th interaction with matter. cm 

n Direction; |n| = 1. 

rii Direction after the i-th interaction with matter. 

v Frequency. Hz 

Vi Frequency after the «-th interaction with matter. Hz 

w Statistical weight. 

Wi Statistical weight after the i-th interaction with matter. 

r/(r, v 1 n) Total emissivity; J^iLo Vi{ r : v > n )- ergs" 1 cm' 3 sr -1 Hz 

r]i(r, v, n) Partial emissivity of light scattered i times. ergs -1 cm -3 sr -1 Hz 

I(r, v, n) Total specific intensity; JDSo -^( r > v i n )- ergs -1 cm -2 sr" 1 Hz 

Ji(r, i/, n) Partial specific intensity of light scattered i times. ergs" 1 cm" 2 sr -1 Hz 

£(r, v, n) Total radiative intensity defined in § 4.4. 2| ; J2ilo £*( r ) v \ n )- ergs" 1 sr" 1 Hz -1 

£j(r, v, n) Partial radiative intensity of light scattered i times. ergs" 1 sr" 1 Hz -1 

I Length. cm 

r(r, v, n; I) Optical depth from r to r + In; J Q dl' x(r + I'n, v, n). 

Too(r, v, n) Optical depth from r to infinity; limj^oo r(r, v, n; /). 

x(r, v, n) Linear extinction coefficient; n + a cm" 1 

<r(r, v, n) Linear scattering coefficient. cm" 1 

«(r, ^, n) Linear absorption coefficient. cm" 1 

a(r, i/, n) Single-scattering albedo; cr/x- 

S(r; i/', n'; v, n) Scattering outcome function, defined by equation [j]. sr" 1 Hz -1 

^(r; v; n', n) Scattering phase function, normalized according to equation |[ sr" 1 



those working in this field. The principal contribution of this paper is to formalize and publish the algorithm, 
so that it can be better understood and so that future researchers can avoid the wasted effort of independent 
rediscovery. 

In §2 we present the restricted problem along with a formal (but computationally intractable) solution. In 
§3 we briefly discuss previous approaches to this problem. In §4 we present the algorithm and a number of 
optimizations and variations. In §5 we discuss important details of our implementations. In §6 we present two 
test cases. In §7 we quantify the efficiency gains due to two important optimizations. In §8 we briefly outline 
how the algorithm might be extended to include time dependence and polarization. In §9 we summarize our 
main results. 

2. THE RESTRICTED PROBLEM 

The restricted problem requires that ?7o(r, ^, n), xir^jn), a(r, v, n), and £(r; z/, n'; l>, n) be specified a 
priori. Table [I] defines our notation, which closely follows that of Mihalas (1978), except for the scattering 
phase function: we use $ instead of g for the phase function, to avoid confusion with the asymmetry parameter 
of the Henyey-Greenstein phase function, and we normalize the phase function according to equation [| so the 
phase function has units of sr" 1 and an isotropic phase function has a value of l/4ir sr" 1 everywhere. Note that 
the emissivities and specific intensities are per unit solid angle and the extinction, scattering, and absorption 
coefficients are per unit length. We denote by Xi the quantity x when only light scattered i times is considered. 
For example, rjo is the emissivity of unscattered light, that is, the true emissivity, and rji is the emissivity 
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of singly-scattered light. For convenience, we also introduce the single-scattering albedo a and the scattering 
outcome function £ The albedo a has its conventional definition as the probability that an interaction with 
matter results in a scattering rather than an absorption. The scattering outcome function £ is defined by 



r](r, v,n) = r) (r, v, n) + / dv' / dd' cr(r, z/,n')£(r;z/,n';z/,n)/(r, z/,n' 

JO J4tt 



) (1) 



That is, £(r; v' , n'; v, n) is the contribution to the emissivity rj(r, v, n) when intensity I(r, v' ', n') is scattered. We 
explicitly separate a and £ to distinguish the probability distribution of scattering events from the probability 
distribution of the outcomes of those events given that they have occurred. 

Scattering does not necessarily conserve energy, but it conserves photons. Considering the specific intensity 
as being carried by photons of energy hv, we can identify £(r; v' , n'; v, n)(z/ jv) as the joint probability distri- 
bution that upon scattering a photon with frequency v' and direction n' becomes a photon with frequency v 
and direction n. Thus, integrating over all initial or final states, we have 

/ di//dft'£(r;z/,n';z/,n)(z//z/) = / dv I dfL £(r; z/, n'; v, n)(i//v) = 1, (2) 

JO J in Jo J iir 

for any v and n or n' and n'. As an example of £, consider coherent scattering (y' = v) with a phase function 
$(r; v\ n', n) normalized as 



/ dVt' $(r; v\ n', n) = j dQ $(r; u; n', n) = 1; 



(3) 



we then have 

£(r; u', n'; v, n)(i//i/) = - i/)$(r; u; n', n). (4) 

The standard formal solution of the equation of radiation transfer gives the intensity Ii(v, v, n) of photons 
that have undergone i scatterings as 

J<(r, i/, n) = /" (fl % (r + /n, i/, n) e - T ( r '^ n; '), (5) 

J l— — OO 

where 

r(r, v, n;Z)= dl' x( r + l' n , l/ ,n). (6) 
Jo 

To evaluate the first integral, we need the emissivities r/i(r, ^, n). The true emissivity 770 (r, v, n) of unscattered 
light is specified a priori. The emissivity ry,(r, u, n) of light that has undergone i scatterings (where i > 0) is 
given by 



%(r,i/,n)= / di/ / df2' cr(r, j/', n')£(r; 1/, n'; 1/, n)7j_i(r, i/', n'). 

JO J47T 



(7) 



In general, equations ^| and ^ form an infinite sequence of coupled integral equations. Approximate solutions 
can be obtained by directly evaluating a few terms and ignoring the rest. The single-scattering-plus-attenuation 
approximation consists of directly evaluating Ji, the intensity of single-scattered light, and ignoring multiply- 
scattered light (I2, -Z3, etc.); this is often adequate for optically-thin problems. However, in optically-thick 
problems, multiply-scattered light can be important, yet direct calculation of the intensity Ii for a single set of 
values of r, u, and n requires the evaluation of a 3i + 1 dimensional integral, and rapidly becomes intractable; 
the Monte Carlo algorithm addresses this intractability. 
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3. PREVIOUS WORK 



The restricted scattering problem described in the preceding section has been tackled many times using 
semi- analytic and Monte Carlo methods. 

Semi-analytic methods include the discrete-ordinate method (Chandrasekhar 1960; Flannery, Roberge, & 
Rybicki 1980; Stamnes et al. 1988) and "doubling" and "adding" methods (van de Hulst 1963; Hansen & Travis 
1974; de Haan, Bosma, & Hovenier 1987). They are well-suited to problems with plane-parallel or spherical 
geometries, but are very difficult to extend to arbitrary geometries. Additionally, these methods become less 
efficient when the scattering phase function is sharply peaked (Escalante 1994), as appears to be the case for 
interstellar grains in the visible and ultraviolet. 

Monte Carlo methods, in which photon trajectories are simulated probabilistically, offer more flexibility. 
Witt (1977), Hillier (1991), Whitney (1991), Fischer, Henning, & Yorke (1994), Watson (1994), Code & Whitney 
(1995), Henney & Axon (1995), and Bjorkman (1997) have presented descriptions of various Monte Carlo 
methods, and Cashwcll & Everett (1959) have described general techniques that can greatly increase the 
efficiency of Monte Carlo methods. Monte Carlo methods have been applied to astrophysical problems many 
times over the last several decades; the earliest application we can find is an investigation of the transfer 
of resonance-line radiation in plane-parallel slabs by Avery & House (1968). Faster computers have allowed 
increasingly complicated geometries to be explored in recent years. 



The algorithm consists of two logically separate parts. Part I follows the course of many pseudo-photons 
as they scatter through the system, producing a sample of "interaction events". These interactions events 
are characterized by the quantities r^, t^, and (the location of the i-th interaction with matter and the 
frequency and direction after the i-th interaction with matter) and by a statistical weight w. We first present a 
naive algorithm for Part I that straightforwardly simulates the underlying physical processes, and then present 
several variations and optimizations, which trade computational efficiency for conceptual complexity. Part II 
makes Monte Carlo estimates of derived quantities, such as the emergent and mean intensities. It is based on 
the observation that the weighted values of the interaction events generated in Part I form unbiased samples 
of the photon emissivities rji/hv and the combinations li\jhv. This algorithm extends those we have seen 
published in its use of forced scatterings, forced interactions, and, in particular, its approach to estimating the 
emergent intensity ( "forced escapes" ) . 



We begin by considering the normalization of the problem. The general radiation transfer problem is non- 
linear, as it includes coupling of radiation and matter. However, the restricted problem is linear, and we can 
impose a convenient normalization on 770 and all derived quantities. Again, considering emissivity in terms of 
photons of energy hv, we can identify r]o(r, v, n)/hu as being proportional to the joint probability distribution 
that a photon is emitted from position r with frequency v and direction n. We can thus impose a normalization 
such that 



so that r]o(r, v, r\)/hv is identically the joint probability distribution of photon emission. All derived quantities 
are thus implicitly normalized by the number of emitted photons. 



4. THE ALGORITHM 



4.1. Normalization 




(8) 



4.2. Part I: Generating The Sample of Interaction Events (Naive Algorithm) 



The sample of interaction events is constructed by simulating the emission and transfer of an adequately 
large number of pseudo-photons; the values of w, rj, v i: and for each pseudo-photon form the sample. The 
algorithm is followed for each pseudo-photon until the pseudo-photon escapes or is absorbed. We first describe 
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a naive version which is a direct analog of the transfer of a real photon (considered as a particle) through the 
real system. 

1: Choose the initial position ro, frequency vq, and direction no from the joint probability distribution 
1o(ro,fo,no)/H. 

2: Initialize the weight: w <— 1. 

3: Initialize the scattering index: i <— 0. 

4: Choose an optical depth from an exponential distribution with a mean of 1. 

If Ti is greater than the optical depth to infinity ^(r^, i/j,n,), then the pseudo-photon escapes: stop. 

5: Otherwise, the pseudo-photon interacts with matter. The interaction occurs after a distance U which is 
given by the implicit equation 

h 

dl' x(r 4 + l'n i: v, TOi) = r(rj, 14, n,; k) = r 4 . (9) 
6: r i+ i <- r t + /n*. 

7: Choose a probability itj from a uniform distribution between and 1. 
If Ui > a(rj + i, nj) then the pseudo-photon is absorbed: stop. 

8: Otherwise, the pseudo-photon is scattered. Choose a frequency v%+i and direction n J+ i from the joint 
probability distribution E(r» + i; z/*, v i+ i, n i+ x){vi/v i+ x). 

9: Increment the scattering index: i <— i + 1. 

10: Continue from step 4. 

4.3. Pari /: Generating The Sample of Interaction Events (Variations and Optimizations) 

We now present two variations and two optimizations. The two variations can simplify the algorithm, easing 
its implementation, at a cost of increasing the variance of the results. The two optimizations can drastically 
reduce the variance in the results at the cost of a small increase in complexity. (This is demonstrated and 
discussed in §|t]). The principal conceptual change in these variations and optimizations is that the statistical 
weight w of the pseudo-photon ceases to be constant. The basis of the variations is weight balancing. An event 
E which should be selected from a probability distribution p(E) can be selected from a different probability 
distribution p'{E) provided p'(E) is non-zero whenever p{E) is non-zero and the statistical weight of the event 
is multiplied by a factor p{E)/p'{E). This technique is useful as it allows an unwieldy probability distribution 
to be replaced by a more straightforward one. The basis of the optimizations is weight splitting. An event E 
with statistical weight w can be replaced by two or more events Ej with statistical weights Wj = wp(Ej\E) 
provided that the Ej fully cover the possible outcomes, that is ~}2p{Ej\E) = 1. This technique can be used to 
reduce the variance in the sample of events. 

4.3.1. Variation: Choosing ro, vq, and no 

In step 1 the initial values of ro, v$, and no are selected from a probability distribution that can be quite 
unwieldy. As an alternative, the values can be selected from different (more easily manageable) distributions 
and then the pseudo-photon can be weighted appropriately. For example, step 1 can be replaced by: 

1': Choose an initial position ro from a uniform distribution whose domain includes all positions r for which 
770 (r, n) is non-zero for some frequency v and direction n. 

Choose an initial frequency vq from a uniform distribution whose domain includes all frequencies v for 
which 770 (r, v, n) is non-zero for some position r and direction n. 
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Choose an initial direction no from a uniform distribution whose domain includes all directions n for 
which r/Q (r, v, n) is non-zero for some position r and frequency v. 

Initialize the statistical weight: w <— 770 ( r o 5 vq, no) /hisy. 

Obviously, variations on this are possible if, for example, some variables can be selected easily from a distri- 
bution but others cannot. All that is required is that the product of the joint probability density for selection 
and the initial statistical weight w be equal to the true joint probability density for emission r/o (i"o , vq, iio)/hvo. 

4.3.2. Variation: Choosing Vi+i and nj+i 

Again, in step 8 the values of Vi+i and n<+i are selected from a probability distribution that can be quite 
unwieldy. As an alternative, the values can be selected from different (more easily manageable) distributionsq 
and then the pseudo-photon can be weighted appropriately. For example, step 8 can be replaced by: 

8': Choose a frequency z^ +1 from a uniform distribution whose domain includes all frequencies v for which 
S(r i+1 ; i/j, n^; v, n) is non-zero for some direction n. 

Choose a direction n^+i from a uniform distribution whose domain includes all directions n for which 
E(i"i + i; i>i, n^; is, n) is non-zero for some frequency v. 

Adjust the statistical weight: w <— u>£(ri+i; Vi, n^; Vi+i, n,+i)(^/^i+i). 

Again, other variations are possible, provided once more that the product of the joint probability density 
for selection and the factor modifying the statistical weight w is equal to the true joint probability distribution 
for scattering S(r i+1 ; v it n^; v i+1 , T^+x)(vi/v i+1 ). 

4.3.3. Optimization: Forced Scatterings 

If forced scatterings are required, step 7 is replaced by the following if i is less than some tuneable parameter: 

7': Split the pseudo-photon into one pseudo-photon that scatters and another that is absorbed. Each pseudo- 
photon has the same set of values of Tj, i/j, and n 3 for j < i and the same r»+i. 

For the pseudo-photon that is absorbed: set the statistical weight to w <— (1 — a(rj+i, Vi, ni))w and stop. 
For the pseudo-photon that is scattered: set the statistical weight to w <— 0(1^+1, i/i,Hi)w and continue. 

The effect of this optimization is that the ramifications of both scattering and absorption are fully explored; 
the pseudo-photon makes appropriately weighted contributions to both the absorbed and scattered events. 
This is particularly important when the albedo is low (see §0). 

4.3.4. Optimization: Forced Interactions 

If forced interactions are desired, step 4 of the algorithm is replaced by the following if i is less than a 
tuneable parameter: 

4': Define the escape probability of the photon along its current trajectory as Pi = e ^ T oc(r i ^i,n i )^ 

Choose a deviate r, from an exponential distribution truncated at ^(r^, z/j, nj) (that is, p(ri) — eT Ti j (1 — 
/%) for < n < Too(rj, z/j, nj) and p(n) = otherwise). 

Split the pseudo-photon into one pseudo-photon that escapes and another that interacts with matter. 
Each pseudo-photon has the same set of values of r 3 , Vj, and n 3 for j < i. 

For the pseudo-photon that escapes: set the statistical weight tow<- ftiW and stop. 

For the pseudo-photon that interacts: set the statistical weight to w <— (1 — Pi)w and continue. 

The effect of this optimization is that the ramifications of both escapes and further interactions are fully 
explored; the pseudo-photon makes appropriately weighted contributions to both the escaping and interacting 
events. This is particularly important when the system is optically thin (see §0). 
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4.4. Part II: Estimates of Derived Quantities 

The sample of interaction events can be used to estimate physically interesting quantities, such as the 
emergent and mean intensity. It is clear from the similarity between the algorithm for Part I and the underlying 
physical process that (a) the weighted values of rj, Vi, and are drawn from the joint probability distributions 
of a photon being emitted or scattered at position with frequency Vi into direction after i scatters and 
(b) the weighted values of r^+i, z/j, and rij are drawn from the joint probability distributions of a photon 
being absorbed or scattered at position r^+i after i previous scatters whilst having frequency V{ and traveling 
in direction n^. These probability distributions are just r]i(r,i',n)/hi' and Jj(r, v, n)x(r, v, x\)/hv) thus the 
weighted values of r^, 2/3, and form unbiased samples of f?i(r, v,n)/hv and the weighted values of r^+i, 2/3, 
and rii form unbiased samples of Jj(r, v, n)x(r, v,r\)/hv. 

We can use standard Monte Carlo techniques on these unbiased samples. If we have a function /(r, v, n) 
defined over a volume V , a frequency interval N , and a solid angle f2, then we have 

W~ x V wf(v h v h m) « [dvfdvfdn Mr.^l/H/fr,^) (10) 



V(r i ,f i ,n i )e(V,JV,n) 



and 



W 1 V w/(r i+ i,i/i,ni) « dV dv dn[Ii(r,v,n)x(r,v,n)/hu]f(r,v,n), 
wta at nl ./v Jn Jn 



(11) 



V(r, + i,i/,,n i )e(V,Ar,f2) 

or equivalently 



W- 1 u;/(ri, !/<_!, !!<_!)« fdV fdu fdn[I^ 1 (r^,n)x(r^,n)/hi^]f(r,iy 1 n). (12) 

V(r < ,i/ i _i,n 4 _i)e(V,iV,n) y w n 

Here, = J^u; is the total statistical weight of all pseudo-photons. The approximations are in the sense 
that the integrals are the means of the sums and, by the central limit theorem, are their limiting values as 
the number of pseudo-photons tends to infinity. Thus, the sums can be used to estimate the integrals. As 
examples, we derive below expressions for the emergent intensity and the mean intensity; other quantities such 
as the heating rate and radiation pressure can be derived similarly. 



4.4.1. Emergent Specific Intensity 

We define a volume V by the translation of a surface S in a direction n from +oo to — oo. We define 
Jj(S, v, n) as the average emergent specific intensity of light scattered i times in direction n from volume V. It 
is given by 

7 i (S,^,ii) = (S-n)- 1 [dV ^(r, i/ > n)e- T -<*'"' n >. (13) 
Jv 

For the scattered emergent intensity, Ii(S, u, n) with i > 0, we can substitute equation ^ for rji to give 

Ji(S,z/,n) = (S-n)- 1 fdvf di/ f dSl' J<_i(r,i/, n>(r,z/, n')£(r; v' , n'; v, n ) e -^(r,»,n) ^ 

JV JO JAit 

We can now use equation [l^ (with N covering all frequencies and Q covering in steradians) to estimate 
Ii(S,v,n) by 

Ii(S,v, n) w (S • n) -1 ^" 1 to hfi-i a(r i: tii_i)E(r 4 ; n^i; v, n)e _T ° o(r " l/ ' n) , (15) 

V(r i ,i/ ( ,n i )e(V r 1 iV,n) 

where, because A and O are all-inclusive, the condition on the sum simplifies to Vr^ € V". 

In some Monte Carlo algorithms, the emergent intensity is estimated by binning the photons that escape 
in step 4 of the algorithm in Part I. Our method of estimating the emergent intensity is very different, in that 
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it is based on the set of interaction events generated in steps 5 to 8. For spherically symmetric problems, the 
two approaches are similar in efficiency. However, for other problems our approach is more efficient by a factor 
of roughly 4tt divided by the size of the bin in solid angle used in the other method; thus, our approach can 
easily be an order of magnitude more efficient for axisymmetric problems and two orders of magnitude more 
efficient for asymmetric problems, provided the viewing direction is fixed. By analogy with forced scatterings 
and forced interactions, we refer to this approach as "forced escapes" . 

We must use a different approach for the unscattered emergent intensity I , as equation ^ is valid only 
for i > and neither equation [l0| nor equation [ll] will give Iq in general. One simple approach is performing 
a logically separate Monte Carlo estimation. That is, selecting values of r' from a uniform distribution that 
includes all positions for which 770 (r, v, n) is non-zero for some v and n and estimating Iq by 

J (S, v, n « S • n)- 1 r ° £V — . (16) 

V L%(ro,^,n)/^ J 

Other Monte Carlo schemes are possible and can take advantage of knowledge of the properties of 770 . The two 
logically separate Monte Carlo estimations can be combined under certain circumstances. 



4.4.2. Emergent Radiative Intensity 

When the system is unresolved, for whatever reason, the specific intensity can no longer be measured. In 
this case, a more useful quantity is the radiative intensity C, which in more familiar terms is the luminosity per 
unit solid angle L = dL/dft, where L is the total luminosity of the system. This quantity is defined and used 
more often in physics than astronomy, although even in physics is it not common. The radiative intensity is 
related to the observed physical flux F by C — Fd 2 , where d is the distance to the source, and to the specific 
intensity I by C — f s IdS = SI(S), where the surface S includes the whole system. Like the specific intensity, 
the radiative intensity has the convenient property of being independent of distance. We can estimate C using 
the same equations as for the specific intensity, but letting the volume of integration extend over the whole 
system and omitting the division by the area S • n. 



4.4.3. Mean Intensity 

The mean value Ji(V,N) of the mean intensity in a volume V and frequency range N of light scattered i 
times is 

Ji(Y,N) = ^ [dV'fdv' Ji(r',v') (17) 



VN 
1 

AkVN 



hnlfh'lf '<<'»'>■ < 18 > 

We can use equation |ll| to estimate Ji by 



1 v—* whvi 

w^™4*vn £ xfr^Ur (19) 

VO,,r 1+1 )e(Af,V) AK 11 11 

Note that while Ii is obtained at a single frequency v, Ji is obtained in a frequency range N. This method 
of calculating the mean intensity seems to be widely used by other Monte Carlo methods, and this particular 
formulation does not enjoy any advantage of efficiency. 



4.5. Estimation of Uncertainty 

In addition to calculating the value of a derived quantity from the whole set of pseudo-photons, we can 
partition the set into M equal sub-sets and obtain M independent values. The uncertainty in the value 
calculated from the whole set can be estimated as the standard deviation of the M values from the M sub-sets 
divided by y/M. 
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5. IMPLEMENTATION 
5.1. Specifying a Problem 

One of the advantages of Monte Carlo methods is their generality and flexibility. For example, to specify 
a problem, our implementations require only five subroutines: to emit a pseudo-photon, to determine the 
extinction coefficient, to determine the albedo, to determine the location of the sentinel surfaces (see below), 
and to scatter a pseudo-photon. 

5.2. Unifying Part I and Part II of the Algorithm 

The algorithm as described above generates a sample of interaction events in Part I and then derives 
quantities from weighted sums of functions of these interaction events in Part II. Since a simulation can 
require many millions of pseudo-photons, a naive implementation would require the storage of many millions 
of interaction events. Our implementations unify Parts I and II, forming the sums as the interaction events are 
generated. This complicates the implementations, but relieves them of the need to store the entire sample of 
interaction events. 

5.3. The Integrator for Optical Depth 

There are many ways to solve equation ^| for the distance li to the optical depth r^, but we have found it 
convenient and efficient to treat it as a differential equation and integrate. That is, we calculate increments 
Al and At to I and r and iterate until we satisfy the boundary condition r = rj, in which case we have 
solved for I = ij. We use a fourth-order Runga-Kutta integrator, which in this case reduces to Simpson's 
rule. Our integrator uses step doubling to estimate the fractional and absolute errors and adapts the step 
size appropriately. If a step fails, it is tried again using a mid-point integrator; this requires no further 
integrand evaluations and so is cheap; furthermore, it allows the integrator to integrate up to and away from 
discontinuities. We have found that this integrator is faster than higher-order Runga-Kutta or Gaussian 
integrators for the extinction distributions we have encountered. We believe that this is because we typically 
require fractional precisions of only 10 or less. The situation is reversed if much higher precisions are imposed. 

We integrate in a piece-wise manner, breaking the integration at "sentinel surfaces" which correspond to 
discontinuities and sharp features such as the equatorial plane of a thin disk. Thus, we are typically integrating 
until we satisfy one of two boundary conditions: until either r = tj or I = L, where L is the distance to the 
next sentinel surface. If the first condition is satisfied, we stop the integration; if the second is satisfied, we 
calculate the distance to the next sentinel surface and continue. This approach allows the integrator to handle 
discontinuities gracefully and ensures that the step size does not become so large that the integrator steps over 
sharp features without noticing them. 

When the integrator is close to one of the boundary conditions, it will often over-shoot. If I over-shoots, we 
adjust the step size so that Al = L — I. If r over-shoots, we adjust the step size so that Al = (r — t.;)/x where 
X is the extinction at the mid-point of the interval. For the initial step we set Al = min(r, T c h)/x where \ is 
extinction at the initial point and r c h is a characteristic optical depth in the problem (often 1); if \ is zero, we 
chose Al = L — I. 

We also use this integrator to determine t^. We do this simply by replacing the boundary condition on 
r with one that / must reach a bounding sphere enclosing the volume in which \ 7^ 0. If we are not forcing 
interactions, we combine steps 4 and 5 of the algorithm into one integration by integrating until either r = Tj 
or I reaches the bounding sphere. 

5.4. Parallelism and Estimation of Uncertainty 

Monte Carlo algorithms are embarrassingly parallel and require very little inter-process communication. 
Once they have started, individual processes only need to communicate at the end of the run to average their 
estimations. This property makes Monte Carlo codes ideal for loosely-coupled clusters of general-purpose com- 
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puters connected by ordinary networks. We have implemented parallelism using the MPICH implementation 
(Gropp et al. 1996) of the Message Passing Interface (MPI). 

The only significant problem is generating different sequences of pseudo-random numbers in each process. 
We use distinct instances of the parallel linear-congruential generators described by L'Ecuyer & Andres (1996). 
Each instance generates distinct sequences with periods of 2 . 

As described above, partitioning of the set of interaction events provides a simple means to empirically 
gauge the errors in the combined estimations; we have used the natural partition that occurs when running 
in parallel to implement this. Even when using only a single-processor or dual- processor computer we often 
employ of order 10 parallel processes and thereby obtain an estimate of the error in the final results; the 
overhead is minimal. 



5.5. The Henyey-Greenstein Phase Function 

The phase function most commonly used for dust scattering problems (when polarization is ignored) is the 
Henyey-Greenstein phase function (Henyey & Greenstein 1941), which is given by 

$hg(^; n, n') = (4tt)- 1 (1 - g 2 v )(l + g 2 . - 2g^ s )-^ 2 , (20) 

where ji s = n • n' is the cosine of the scattering angle 9 S and the asymmetry parameter g v is the mean value of 
/i s . With this choice of scattering function, step 8 of the algorithm is considerably simplified, as /i s and hence 
the scattering angle 9 S can be obtained from 

J 1 - 2a for g v = 

MS \ (2<?,)- 1 (l+^-(l-^) 2 (l+.9,(l-2a))- 2 ) for.g^0 ' 1 > 

where a is a deviate drawn from a uniform distribution between and 1. The azimuthal scattering angle 4>s is 
uniformly distributed between and 27r. 



6. TEST CASES 



Scattering algorithms seem to be tricky to implement correctly, and simple but non-trivial test cases are 
difficult to find in the literature. For this reason, we present two simple test cases. These test cases are 
interesting in that they are non-trivial, yet the first few intensities can be obtained by direct numerical in- 
tegration. The first test case tests the normalization of a point source and the ratio of the unscattered and 
singly-scattered intensities. The second test case tests the normalization of a pencil beam of light and the 
ratios of the unscattered, singly-scattered, and doubly-scattered intensities. 



6.1. Point Source and Slab 



Consider a mono-chromatic, isotropic point source illuminating a plane-parallel slab that is infinite in the 
x and y directions but has a finite optical depth T in the z direction. Scattering is coherent. What is £oo, 
the radiative intensity seen by an observer at infinity, as a function of 6 = arccos/x, the angle from the +z 
direction? (We discuss the radiative intensity in § 4.4. 2\ .) 



To fix the geometry, we can take the true emissivity and extinction distribution to be 

(22) 

\ 0<Z<1 • (23) 
otherwise 

This problem is similar to that considered by Henney (1998) in his §2.1.2, except that here we consider a point 
source rather than an infinite sheet source (and we also normalize $ differently). Therefore, the radiative 
intensities for this test case will simply be his specific intensities multiplied by a factor of and normalized; 
this is a general technique for transforming between finite unresolved and infinite resolved sources. Thus, 
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£00,0 an d £00,1, the unscattered and singly-scattered emergent radiative intensities seen from infinity, are given 
analytically by 

£oo,o(M) = (24) 

au,A(u) f 1 f 27T 1 - p-^-aO/w' 
£oc,i(/x) = /d/x' / # ^(/i,//',^)) ; (25) 

where 

f e ~ T/M if > 

ll II jJL < 

Ms(M,M',0) = Wi' + (l-M 2 ) 1/2 (l-M ,2 ) 1/2 cos^. (27) 

Spatially unresolved quantities seen from infinity are independent of the precise form of the extinction 
coefficient distribution within the slab, provided we conserve the plane-parallel symmetry and the optical 
depth through the slab. Therefore, when solving this problem analytically, we can treat the slab as uniform. 
When solving it with the Monte Carlo algorithm, we can choose an extinction distribution that provides a more 
severe test of the integrator, such as 

f T[l/2 + sin 2 (27rz)] if < z < 1 
[ otherwise 

Figure [j] compares the direct and Monte Carlo estimates of the emergent radiative intensities for a slab with 
T = 2, a = 0.5, and a Henyey-Greenstein phase function with g — 0.5. The values of £00,0 an d £00,1 calculated 
by the two methods agree to 10 -4 or better, which is commensurate with our estimates of the precisions of the 
two calculations. Selected values of £ 00j 0, £00,1, £00,2, an d £00, >2 ar e given in Table |. 

6.2. Pencil Beam and Slab 

For our second test case, we use the same slab but replace the point source with a pencil beam along the 
+z axis. We can take the true emissivity to be 

Tk>(r,n) = 6(T)5(ii). (29) 

Equations || and fj] can be used to show that £00,0, £oo,i 5 and £00,2, the unscattered, singly-scattered, and 
doubly-scattered emergent radiative intensities seen from infinity, are given analytically by 

£oo,o(M) - e- T B(fi) (30) 

£oo,i(M) = aC(M)$(M) T ^L(l-e- r ( 1 -^) ) (31) 

r T r+ l r 2 ' !T ,, ... 1 _ p -t\t,ii'){i-ii')/\ij.'\ 

£oo,2(M) = a 2 dr dfi' / # $(^y,<«)$(AOe- r e~ r (r '~ Al)/H — ; ; , (32) 

Jo J-i Jo 1 — M 



where 
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8 

Fig. 1. Results for the point source and slab test case described in § |6.l[ The figure shows the emergent radiative 
intensity Cao(9). The lines are calculated directly, with the solid line showing unscattered light £00,0 and the dashed 
line showing singly-scattered light £00,2- The symbols are calculated using the Monte Carlo algorithm, with circles o 
showing the radiative intensity of unscattered light £00,0, squares □ showing the radiative intensity of singly-scattered 
light £00,1, triangles A showing the radiative intensity of double-scattered light £00,2, and diamonds o showing the 
radiative intensity of more-than-doubly-scattered light £00, >2 = 
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Fig. 2. Results for the pencil beam and slab case described in 

Coo{0). The lines are calculated directly, with the dashed line showing singly-scattered light £oo,i 
line showing doubly-scattered light £oo,2- Fhe symbols are calculated using the Monte Carlo algorithm, with circles o 
showing the radiative intensity of unscattered light £00,1, squares □ showing the radiative intensity of singly-scattered 
light £oo,2, triangles A showing the radiative intensity of doubly-scattered light £00,2, and diamonds o showing the 
radiative intensity of more-than-doubly-scattered light C 
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Fig. 3. Contour plots of the relative errors in three trials of the slab geometry as functions of the number of forced 
scatterings and forced interactions. Trial (a) has T — 2 and a — 0.5, trial (b) has T = 0.1 and a = 0.5, and trial (c) has 
T — 2.0 and a = 0.1; in all cases g = 0.5 and 6 — 45°. Contours are spaced by factors of 2. The hatched regions have 
the lowest errors. 



These results are again most easily derived using the transformation mentioned in the previous section. Note 
that as the singly-scattered light can be written in a closed form, the doubly-scattered light is given by a single 
3-dimensional integral. 

Figure || compares the direct and Monte Carlo estimates of the emergent radiative intensities for a slab with 
T = 2, a = 0.5, and a Henyey-Greenstein phase function with g = 0.5. The values of £oo,i and £00,2 calculated 
by the two methods agree to 5 x 10 -4 or better, which is commensurate with our estimates of the precisions of 
the two calculations. Selected values of £00,0, A», 1, £00,2, and £00, >2 are given in Table ||. 



7. THE IMPORTANCE OF FORCED SCATTERINGS AND INTERACTIONS 



To illustrate the importance of forced scatterings (§4.3.3) and forced interactions (§4.3.4), we have examined 



the relative errors in the emergent intensity in simulations of the slab geometry (§|6|) as functions of the number 
of forced scatterings and the number of forced interactions. In order to make the comparison fair, we ran 
each calculation for the same amount of processor time. The relative errors were estimated by partitioning the 
sample of pseudo-photons into 48 sub-samples. 

We ran three trials with different values of T and a. Trial (a) has T — 2.0 and a — 0.5, the same as the 
test case considered in §^, and is moderately optically-thick, with singly-scattered and doubly-scattered light 
dominating more-than-doubly-scattered light. Trial (b) has T = 0.1 and a = 0.5, and is optically-thin, with 
singly-scattered light dominating multiply-scattered light. Trial (c) has T — 2.0 and a = 0.1, and is moderately 
optically-thick, but has a low albedo, so singly-scattered light also dominates multiply-scattered light. In all 
cases we kept 8 = 45° and g = 0.5. 

Figure || shows contour plots of the relative errors for these trials. The best results for trial (a), with T = 2 
and a — 0.5, are obtained in a broad region with at least 4 forced scatters and 3 forced interactions. The best 
results for trial (b), with T = 0.1 and a = 0.5, were obtained for at least 2 forced interactions and at least 2 
forced scatters. The best results for trial (c), with T = 2.0 and a — 0.1, were obtained for at least 2 forced 
interactions and at least 3 forced scatters. In all cases, however, the optimal region is roughly L-shaped. That 
is, excessive numbers of either forced scatters or forced interactions have little effect, but excessive numbers of 
both start to increase the error again. 

The best errors for trials (a), (b), and (c) are roughly 7, 90, and 30 times smaller than for a naive imple- 
mentation with no forced scatters or forced interactions Since the relative error in Monte Carlo calculations 
decreases as the square root of the computational effort, these optimizations represent savings in time of factors 
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of roughly 50, 8000, and 1000. 

These results can be understood qualitatively as follows. On average, the effect of forcing scatterings and 
interactions is to keep the pseudo-photon in the system for longer, which has two effects: (i) the pseudo-photon 
is more likely to contribute to derived quantities that depend on scattered light, which reduces the variance in 
these quantities, and (ii) the computational cost per pseudo-photons increases, so fewer can be followed, which 
increases the variance in derived quantities. In the cases examined here, the former is more important to begin 
with; for moderate numbers of forced scatterings and interactions, making each pseudo-photon more useful 
compensates for the smaller total number of pseudo-photons, and the variance is reduced. This is especially 
important if the system is optically thin (trial b) or if the albedo is low (trial c), as pseudo-photons can easily 
leave the system or be absorbed, and thereby make relatively little contribution to derived quantities. 

However, if the number of both forced scatterings and interactions is excessive, each pseudo-photon is forced 
to remain in the system for many scatterings. This may reduced the variance of quantities derived from very 
highly scattered light, but it reduces the total number of pseudo-photons that can be considered, and thereby 
increases the variance of derived quantities that are dominated by moderately scattered light. Furthermore, the 
optimal region is L-shaped because, for the optical depths and albedos considered here, both excessive numbers 
of forced scatterings and excessive numbers of forced interactions are required to keep each pseudo-photons in 
the system for many scatterings; if one or the other is not excessive, the pseudo-photon is absorbed or escapes. 

It is impossible to give a universally applicable rule to select a priori the optimal number of forced scat- 
terings and interactions, as these quantities depend on the geometry and also on the particular quantity being 
calculated. However, the discussion in the previous paragraph implies that the optimal values for both will 
probably be roughly equal. Our experience suggests that two to four forced scatterings and interactions might 
be a useful rule of thumb; this number is not too far from optimal for any of the cases we have investigated. 

Naive Monte Carlo algorithms are often considered ill-suited to optically-thin problems, as the overwhelming- 
majority of photons escape without interacting. However, is not the case for more sophisticated algorithms 
that force the first few interactions, as our optically-thin trial demonstrates. A specialized single-scattering 
plus attenuation calculation would probably calculate the emergent singly-scattered intensity more efficiently, 
but a Monte Carlo calculation may well be more flexible and can directly account for multiply-scattered light. 

8. TIME-DEPENDENCE AND POLARIZATION 

The algorithm as presented covers the time-independent transfer of unpolarized light. Both of these restric- 
tions can be lifted quite easily. 

The simplest way of including general time-dependence is to sample the range of times of emission, keep 
track of the travel time of the pseudo-photon within the system, and account for the travel time when computing 
quantities of interest. If only the source varies, it is probably more efficient to construct a "response function" 
for each quantity of interest which can be convolved with changes in the brightness of the source to predict the 
time- variation of that quantity. 

Perhaps the simplest way to include polarization is to generalize the scalar statistical weight w to a vector 
statistical weight w = {w[ 1 wq,wjj,wv) 1 whose components correspond to the four components of the Stokes 
vector. The components will in general not be conserved upon scattering, but will transform according to the 
adopted scattering matrix. See, for example, Hillicr (1991) and Whitney (1991). 

9. SUMMARY 

We have presented and discussed a Monte Carlo algorithm for a restricted class of scattering problems. 
Our main contributions have been to formalize this algorithm, to present two test cases, and to investigate its 
efficiency Important features of the algorithm are forced escapes, forced scatterings, and forced interactions. 
Using forced escapes to estimate the emergent intensity can be orders of magnitude more efficient than naively 
binning escaping pseudo-photons. Also, using a small number (two to four) of forced scatterings and forced 
interactions can significantly improve the efficiency of the algorithm for many problems, overwhelmingly so for 
optically-thin problems or those where the albedo is small. 

We are grateful to Karl Stapelfcldt and John Krist for serving as crash-test dummies for several years, to 
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Jon Bjorkman, Kenny Wood, and Barbara Whitney for many useful discussions on Monte Carlo techniques 
and applications, and to an anonymous referee. This work was supported by CONACyT project 27570E. 
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